/**
 * @file Multigrid_Solver_1D.h
 * @brief A simple implementation for Multigrid Solver; 
 * It deals with 1-D Possion Equation on [0,1] with Dirichlet Boundary Condition;
 * @author Linhuo (1066117119@qq.com)
 * @version 1.0
 * @date 2021-06-07
 * 
 * @copyright Copyright (c) 2021  Linhuo
 * 
 */
#ifndef _PRO2_MULTIGRID_SOLVER_1D_H_
#define _PRO2_MULTIGRID_SOLVER_1D_H_
#include "Multigrid_Solver.h"

/**
 * @brief A simple implementation for Multigrid Solver; 
 * It deals with 1-D Possion Equation on [0,1] with Dirichlet Boundary Condition;
 */
class MG_1D : public Multigrid_Solver
{
public:
    void Initialize_u();
    void Initialize_BC(Function _BCF);
    void InitializeRhs(Function _F);
    VectorXd SmoothingforVC(VectorXd _u, VectorXd _f, double _h);
    VectorXd ResidualSolver(const VectorXd _u, const VectorXd _f, const double _h);
    VectorXd RestrictSolver(VectorXd _r);
    VectorXd IntpSolver(VectorXd _r);
    VectorXd PrepareFMG(VectorXd _f);
    int get_now_num_grid(VectorXd _f) { return (_f.rows() - 1); };
};

VectorXd MG_1D::PrepareFMG(VectorXd _r)
{
    VectorXd __u;
    int N = get_now_num_grid(_r);
    __u.setZero(N + 1);
    __u(0) = _BC(0, 0);
    __u(N) = _BC(1, 0);
    return __u;
};

void MG_1D::Initialize_u()
{
    int N = _num_of_grid;
    _u_origin.setZero(1 + N);
    _u_origin(0) = _BC(0, 0);
    _u_origin(N) = _BC(1, 0);
};

void MG_1D::Initialize_BC(Function _BCF)
{
    double a[1] = {0};
    _BC.setZero(2, 1);
    _BC(0, 0) = _BCF(a);
    a[0] = 1;
    _BC(1, 0) = _BCF(a);
};

void MG_1D::InitializeRhs(Function _F)
{
    int num = _num_of_grid;
    double a[1] = {0};
    _f_origin.resize(num + 1);
    for (int i = 1; i < num; i++)
    {
        a[0] = 1.0 * i / num;
        _f_origin(i) = _F(a);
    };
    _f_origin(0) = _BC(0, 0);
    _f_origin(num) = _BC(1, 0);
};

VectorXd MG_1D::SmoothingforVC(VectorXd _u, VectorXd _f, double _h)
{
    VectorXd _u_temp;
    int N = _u.rows() - 1;
    _u_temp.setZero(N + 1);
    /// Since the first and the last term of _u is ghost cell，we keep them the same in _u_temp.
    _u_temp(0) = _u(0);
    _u_temp(N) = _u(N);
    double _t = 0;
    double omega = 2.0 / 3;
    switch (SmoothType)
    {
    case GaussSeidel:
        for (int i = 1; i <= N - 1; i++)
        {
            _u_temp(i) = (_u_temp(i - 1) + _u(i + 1) - _h * _h * _f(i)) / 2.0;
        };
        break;
    case WeightedJacobi:
        for (int i = 1; i <= N - 1; i++)
        {
            _t = (_u(i - 1) + _u(i + 1) - _h * _h * _f(i)) / 2.0;
            _u_temp(i) = (1.0 - omega) * _u(i) + omega * _t;
        };
        break;
    default:
        break;
    };
    return _u_temp;
};

VectorXd MG_1D::ResidualSolver(const VectorXd _u, const VectorXd _f, const double _h)
{
    VectorXd _u_temp;
    int N = _u.rows() - 1;
    _u_temp.setZero(N + 1);
    for (int i = 1; i <= N - 1; i++)
        _u_temp(i) = _f(i) - (_u(i - 1) - 2.0 * _u(i) + _u(i + 1)) / (_h * _h);
    return _u_temp;
};

VectorXd MG_1D::RestrictSolver(VectorXd _r)
{
    VectorXd _u_temp;
    int N = (_r.rows() - 1) / 2;
    _u_temp.setZero(N + 1);
    switch (RestricType)
    {
    case FullWeighting:
        for (int i = 1; i <= N - 1; i++)
        {
            _u_temp(i) = (_r(2 * i - 1) + 2.0 * _r(2 * i) + _r(2 * i + 1)) / 4.0;
        };
        break;
    case Injection:
        for (int i = 1; i <= N - 1; i++)
        {
            _u_temp(i) = _r(2 * i);
        };
        break;
    default:
        break;
    };
    _u_temp(0) = _r(0);
    _u_temp(N) = _r(2 * N);
    return _u_temp;
};

VectorXd MG_1D::IntpSolver(VectorXd _r)
{
    VectorXd _u_temp;
    int N = (_r.rows() - 1) * 2;
    _u_temp.setZero(N + 1);
    switch (IntpType)
    {
    case Linear:
        for (int i = 1; i <= N - 1; i = i + 2)
            _u_temp(i) = (_r((i - 1) / 2) + _r((i + 1) / 2)) / 2.0;
        for (int j = 0; j <= N; j = j + 2)
            _u_temp(j) = _r(j / 2);
        break;
    case Quadratic:
        // _u_temp(1) = _r(0) / 4.0 + 5.0 * _r(1) / 8.0 + _r(2) / 8.0;
        // for (int i = 3; i <= N - 3; i = i + 2)
        // {
        //     _u_temp(i) = 3.0 * _r((i - 1) / 2) / 16.0 + 5.0 * _r((i + 1) / 2) / 8.0 + 3.0 * _r((i + 3) / 2) / 16.0;
        // };
        // _u_temp(N - 1) = _r(N / 2 - 2) / 8.0 + 5.0 * _r(N / 2 - 1) / 8.0 + _r(N / 2) / 4.0;
        // for (int j = 0; j <= N; j = j + 2)
        //     _u_temp(j) = _r(j / 2);
        _u_temp(1) = _r(0) / 2.0 + _r(1) / 2.0;
        for (int i = 3; i <= N - 3; i = i + 2)
        {
            _u_temp(i) = 9.0 * _r((i - 1) / 2) / 16.0 + 9.0 * _r((i + 1) / 2) / 16.0 - 1.0 * _r((i + 3) / 2) / 16.0 - 1.0 * _r((i - 3) / 2) / 16.0;
        };
        _u_temp(N - 1) = _r(N / 2 - 1) / 2.0 + _r(N / 2) / 2.0;
        for (int j = 0; j <= N; j = j + 2)
            _u_temp(j) = _r(j / 2);
        break;
    default:
        break;
    };
    _u_temp(0) = _r(0);
    _u_temp(N) = _r(N / 2);
    return _u_temp;
};

#else
//DO NOTHING
#endif